rm(list = ls())

europeLink="https://github.com/Spatial-Data-Analytics-DACSS-690D/HW1-Interactive-Visualization/raw/refs/heads/main/WORLD/europe_8857.gpkg"
library(sf)
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
st_layers(europeLink)

Reading the layers:

europe <- st_read(europeLink, layer = "continent")
## Reading layer `continent' from data source 
##   `https://github.com/Spatial-Data-Analytics-DACSS-690D/HW1-Interactive-Visualization/raw/refs/heads/main/WORLD/europe_8857.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 37 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -2668477 ymin: 3480747 xmax: 3207182 ymax: 7761460
## Projected CRS: WGS 84 / Equal Earth Greenwich
europemob_ddm <- st_read(europeLink, layer = "mobiles_ddm")
## Reading layer `mobiles_ddm' from data source 
##   `https://github.com/Spatial-Data-Analytics-DACSS-690D/HW1-Interactive-Visualization/raw/refs/heads/main/WORLD/europe_8857.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 14480 features and 0 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -2662224 ymin: 3515365 xmax: 3197662 ymax: 7749696
## Projected CRS: WGS 84 / Equal Earth Greenwich
europemob_psm <- st_read(europeLink, layer = "mobiles_psm")
## Reading layer `mobiles_psm' from data source 
##   `https://github.com/Spatial-Data-Analytics-DACSS-690D/HW1-Interactive-Visualization/raw/refs/heads/main/WORLD/europe_8857.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 37 features and 5 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -1290663 ymin: 4371225 xmax: 2918367 ymax: 7326054
## Projected CRS: WGS 84 / Equal Earth Greenwich

Checking CRS

st_crs(europe)$epsg
## [1] 8857

THE BASE

library(ggplot2)
library(scales)

ggplot() + theme_light() +
  #call map
  geom_sf(data = europe) + 
  # adjust datum 
  coord_sf(datum = st_crs(europe)) + 
  # custom X and Y axes
  scale_x_continuous(labels = label_number(scale = 1/1000, suffix = " km")) +
  scale_y_continuous(labels = label_number(scale = 1/1000, suffix = " km")) +
  labs(x = "Easting", y = "Northing") #labels

# DDM - static

plot1_a=ggplot() + theme_light() +
  #base
  geom_sf(data = europe,fill="white",color="grey80") +
  # layer
  geom_sf(data=europemob_ddm,
          alpha=0.1,#transparency
          shape=20,
          size=0.5,
          color="red")+
  coord_sf(datum = st_crs(europe)) + 
  # custom X and Y axes
  scale_x_continuous(labels = label_number(scale = 1/1000, suffix = " km")) +
  scale_y_continuous(labels = label_number(scale = 1/1000, suffix = " km")) +
  labs(x = "Easting", y = "Northing") #labels

plot1_a

# DDM - interactive

library(leaflet)
# needs geographic CRS
ddm_4326 <- st_transform(europemob_ddm, crs = 4326)

leaflet() %>%
  addTiles() %>%  # base map
  addCircleMarkers(data = ddm_4326,
                   color="red",
                   stroke=0,
                   radius = 0.01,
                   opacity = 0.1)

Almost the same:

plot1_b=leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron)  %>%  # names in english
  addCircleMarkers(data = ddm_4326,
                   color="red",
                   stroke=0,
                   radius = 0.01,
                   opacity = 0.1)
plot1_b

PSM - static

plot2_a=ggplot() + theme_light() +
  #base
  geom_sf(data = europe,fill="white",color="grey80") +
  # layer
  geom_sf(data=europemob_psm,
          shape=20,
          aes(size=mobiles/1000,color=factor(mobiles_outlier)))+
  scale_size_area(
    name = "Mobile Users (x1,000)", # Name for the legend
    max_size = 20) +       # This defines the maximum *radius*
  guides(color="none") + 
  
  coord_sf(datum = st_crs(europe)) + 
  # custom X and Y axes
  scale_x_continuous(labels = label_number(scale = 1/1000, suffix = " km")) +
  scale_y_continuous(labels = label_number(scale = 1/1000, suffix = " km")) +
  labs(x = "Easting", y = "Northing") #labels

plot2_a

# PSM - interactive

# CRS for leaflet
psm_4326 <- st_transform(europemob_psm, crs = 4326)


plot2_b=leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron)  %>%  
  addCircleMarkers(data = psm_4326,
                   # fillColor = "red",
                   # stroke = T, #border
                   # color="red", #border color
                   # fillOpacity = 1,
                   label = ~Country,
                   radius = ~size)
plot2_b

CHORO - static

plot3_a=ggplot() + theme_light() +
  #base
  geom_sf(data=europe,
          aes(fill=mobiles_density_FJ5_cat))+
  coord_sf(datum = st_crs(europe)) +
  scale_fill_brewer(direction = 1,palette = 'RdBu') + 
  scale_x_continuous(labels = label_number(scale = 1/1000, suffix = " km")) +
  scale_y_continuous(labels = label_number(scale = 1/1000, suffix = " km")) +
  labs(x = "Easting", y = "Northing",fill="Mobile Users Density")
  
plot3_a

CHORO - interactive

choro_4326 <- st_transform(europe, crs = 4326)
#needed
choro_4326$mobiles_density_FJ5_cat=ordered(choro_4326$mobiles_density_FJ5_cat)
# First, create a color palette for your categorical variable
pal <- colorFactor(palette = "RdBu",
                   domain = choro_4326$mobiles_density_FJ5_cat,
                   ordered = T,
                   reverse = F)

leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron)  %>% 
  addPolygons(
    data = choro_4326,
    fillColor = ~pal(mobiles_density_FJ5_cat),
    fillOpacity = 0.7,
    color = "white",
    weight = 1,
    opacity = 1
    ) %>%
  addLegend(
    pal = pal,
    values = choro_4326$mobiles_density_FJ5_cat,
    )
# Calculate the center of europe for the home button
europe_center <- st_bbox(choro_4326) %>% 
  as.vector() %>% 
  {c(mean(c(.[1], .[3])), mean(c(.[2], .[4])))}  # Calculate center coordinates


plot3_b=leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron)  %>% 
addPolygons(
    data = choro_4326,
    fillColor = ~pal(mobiles_density_FJ5_cat),
    fillOpacity = 0.7,
    color = "white",
    weight = 1,
    opacity = 1,
    label = lapply(
      paste("<strong>", choro_4326$Country, "</strong><br>",
            "Density: ", choro_4326$mobiles_density),
      htmltools::HTML
    )
  ) %>%
  addLegend(
    pal = pal,
    values = choro_4326$mobiles_density_FJ5_cat,
    opacity = 0.7,
    title = "Mobile Users Density",
    position = "bottomright"
  )%>%
  addEasyButton(
    easyButton(
      icon = "fa-home", 
      title = "Reset to europe View",
      onClick = JS(sprintf("function(btn, map){ map.setView([%f, %f], 3); }", 
                          europe_center[2], europe_center[1]))
    )
  )
plot3_b
saveRDS(plot1_a,file = "plot1_a.rds")
saveRDS(plot1_b,file = "plot1_b.rds")
saveRDS(plot2_a,file = "plot2_a.rds")
saveRDS(plot2_b,file = "plot2_b.rds")
saveRDS(plot3_a,file = "plot3_a.rds")
saveRDS(plot3_b,file = "plot3_b.rds")